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EVALUATION OF THE LATTICE-BOLTZMAN EQUATION SOLVER POWERFLOW 

FOR AERODYNAMIC APPLICATIONS 

DAVID P. LOCKARD*, LI-SHI LUO+, AND BART A. SINGER* 

Abstract. A careful comparison of the performance of a commercially available Lattice-Boltzmann Equa- 
tion solver (PowerFLOW) was made with a conventional, block-structured computational fluid-dynamics 
code (CFL3D) for the flow over a two-dimensional NACA-0012 airfoil. The results suggest that the ver- 
sion of PowerFLOW used in the investigation produced solutions with large errors in the computed flow 
field; these errors are attributed to inadequate resolution of the boundary layer for reasons related to grid 
resolution and primitive turbulence modeling. The requirement of square grid cells in the PowerFLOW 
calculations limited the number of points that could be used to span the boundary layer on the wing and 
still keep the computation size small enough to fit on the available computers. Although not discussed in 
detail, disappointing results were also obtained with PowerFLOW for a cavity flow and for the flow around 
a generic helicopter configuration. 

Key words, lattice Boltzmann method, PowerFLOW, 2D flow past NACA-0012 airfoil, aerodynamic 
simulations, turbulence modeling 

Subject classification. Fluid Mechanics 

1. Introduction. In recent years, methods of the lattice-gas automata (LG A) and the lattice Boltz- 
mann equation (LBE) have become an alternative to conventional computational fluid dynamics (CFD) 
methods for various systems (see collections of papers [16], proceedings [7, 6], reviews [4, 14, 30], and mono- 
graphs [38, 15]). The lattice-gas and lattice Boltzmann methods have been particularly successful in dealing 
with fluid flow applications involving complex fluids and complicated boundaries, such as multi-component 
fluids through porous media [8, 31]. Although the promising potential of the LGA and LBE methods as 
viable CFD tools has been demonstrated in a number of cases of laminar flow with low Reynolds number, 
there is no experience to demonstrate that the method can be applied to turbulent flow in the presence of a 
streamlined body at high Reynolds number. In fact, application of the LBE method to turbulent flows at 
high Reynolds number remains as an area of future development [44, 30] . 

The fundamental philosophy of the lattice gas automata and the lattice Boltzmann equation is to con- 
struct simple models based on kinetic theory that preserve the conservation laws and necessary symmetries 
such that the emerging behavior of these models obeys the desired macroscopic equations. The theoretical 
foundation for using simple kinetic models to simulate hydrodynamic systems (or other complex systems) 
is based upon the observation that in nature macroscopic phenomenon such as hydrodynamics is rather 
insensitive to the underlying details in the microscopic, or even mesoscopic, dynamics. In hydrodynamic 
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systems, details of microscopic dynamics can affect the numerical values of the transport coefficients, but not 
the form of the Navier-Stokes equations. Therefore, as long as these models preserve conservation laws and 
necessary symmetries (e.#., Galilean invariance) satisfactorily, they can be used to simulate hydrodynamic 
systems. The primary design criteria for the development of LGA and LBE models are efficiency and the 
preservation of the physical properties of the system. The LGA and LBE models resemble molecular dynam- 
ics methods or the Boltzmann kinetic equation in concept, but do not have their computational complexity. 
However, some artificial physical limitations are typically placed on the system to improve the efficiency of 
the algorithm. 

The kinetic nature of the LGA and LBE methods leads to the following features that distinguish the 
LGA and LBE methods from any other conventional CFD method. First, the convection operator of the 
LGA or LBE models is linear in phase space, similar to that of the Boltzmann kinetic equation but different 
than the Euler or the Navier-Stokes equations. Second, the pressure is obtained through an equation of state, 
as opposed to solving a Poisson equation in the incompressible Navier-Stokes equations. Third, with the 
Navier-Stokes equations, the constitutive relations are input from empirical data, while with the LGA and 
LBE methods, the constitutive relations emerge as a result of proper modeling of interparticle potentials. 
Fourth, unlike the Navier-Stokes or Euler equations, in which macroscopic conservation laws are discretized, 
the LGA and LBE methods utilize a set of discrete particle velocities such that the conserved quantities are 
preserved up to machine accuracy in the calculations. In addition, the LGA and LBE methods have the 
following computational characteristics: 

1. Intrinsic parallelism due to nearest neighbor data communications of the streaming (convection) 
process and purely local calculation of the collision process; 

2. Ability to handle complex boundaries without compromising computational speed. 

The intrinsic parallelism of the LGA and LBE methods is an appealing feature in light of the clear trend 
to use massively parallel computers. The ability to handle complex boundaries is important to practical 
applications such as flow through porous media. The capability to include model interactions is crucial to 
simulations of some physical systems such as multi-phase or multi-component flows. Furthermore, the LGA 
method has the following features, in addition to the above, due to its Boolean nature: 

1. Unconditional stability; 

2. Memory efficiency; 

3. Suitability for special-purpose computers. 

The theoretical advantages of the LGA and LBE methods are still largely untested in the realm of CFD. 
Moreover, precisely what advantages the LGA and LBE methods exhibit over conventional methods of 
solving the Navier-Stokes equations are expected to be problem dependent. The areas where the LGA and 
LBE methods are more suitable depend on the physical nature of the problems. Unfortunately, due to 
the lack of research effort in this area, a sufficient amount of numerical evidence to verify the theoretical 
advantages of the LGA and LBE methods does not exist, especially in the areas where conventional CFD 
methods have proven to be useful, such as in the field of aerodynamics. 

Recently, a commercial CFD software product, PowerFLOW, based on the lattice-gas and the lattice 
Boltzmann method became available [17]. PowerFLOW has been applied to simulate high Reynolds number 
turbulent external flow over an automobile body in ground proximity [3]. However, PowerFLOW has not 
been subjected to a set of thorough validation tests. The focus of this study is to explore the possibility of 
using PowerFLOW in typical aerodynamic applications of CFD simulations of various configurations. In this 
report, we present test results for a 2-D flow past a streamlined airfoil (NACA-0012) [1] with Reynolds num- 
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bers (Re = UL/v , where U is the free-stream velocity, L is the chord length, and v is the kinematic viscosity) 
ranging from 500 to 500 x 10 6 using PowerFLOW and compare the results with those from a conventional 
CFD solver CFL3D [48, 51, 47, 40, 39, 29]. Our test results regarding PowerFLOW are rather negative: in 
its present state, PowerFLOW, release 3.0 is inadequate for aerodynamic applications. Nevertheless, with 
further development, LGA and LBE methods may prove to be useful tools for aerodynamics. 

This report is organized as follows. In the next section, the basic notion of lattice-gas automata and 
the lattice Boltzmann equation is briefly reviewed. The lattice Boltzmann equation is presented as a special 
finite-difference form of the continuous Boltzmann equation, and the boundary conditions and turbulence 
models in the lattice Boltzmann method are discussed. Sections 3 and 4 briefly describe PowerFLOW and 
CFL3D, respectively. Section 5 presents the simulation results of flow past a NACA-0012 airfoil for Reynolds 
numbers between 500 and 500 x 10 6 . Section 6 discusses the results and concludes this report. 

2. Lattice-Gas Automata and Lattice Boltzmann Equation. 

2.1. Lattice Gas Automata. In a series of articles published in the 1980’s [53], Wolfram showed 
that cellular automata, despite their simple construction, have sufficient complexity to accomplish universal 
computing: that is, beginning with a particular initial state, the evolution of some automaton could realize 
any chosen finite algorithm. Based upon kinetic theory and the previous experience of the Ha.rdy-Pomeau- 
de Pazzis (HPP) model [21, 22] that a two-dimensional square lattice does not possess sufficient symmetry for 
hydrodynamics, Frisch, Hasslacher, and Pomeau [19] and Wolfram [54] independently discovered that a simple 
cellular automaton on a two-dimensional triangular lattice can simulate the Navier-Stokes equations. The 
two-dimensional Frisch-Hasslacher-Pomeau (FHP) model was immediately generalized to a three-dimensional 
LGA model [20]. 

The LGA model proposed by Frisch, Hasslacher, and Pomeau, and Wolfram evolves on a two-dimensional 
triangular lattice space. The particles have momenta that allow them to move from one site on the lattice 
to another in discrete time steps. A particular lattice site is occupied by either no particle or one particle 
with a particular momentum pointing to a nearest neighbor site. Therefore, at most a lattice site can be 
simultaneously occupied by six particles, hence this model is called the 6- velocity model or FHP-I model. The 
evolution of the LGA model consists of two steps: collision and advection. The collision process is partially 
described in Fig. 1. For example, two particles colliding with opposite momenta will rotate their momenta 
60° clockwise or counter-clockwise with equal probability. In Fig. 1, we do not list those configurations which 
can be easily obtained by rotational transformation, and which are invariant under the collision process. The 
particle number, the momentum, and the energy are conserved in the collision process locally and exactly. 
(Because the FHP model has only one speed, the energy is no longer an independent variable, it is equivalent 
to the particle number. However, for multi-speed models, the energy is an independent variable.) 

The evolution equation of the LGA can be written as: 

(2.1) 7l(x{x &<xSt) t *T St) — Tl^(«Z?, t) C a, i 

where n a is the Boolean particle number with the velocity e a , C a is the collision operator, x is a vector in 
the lattice space with lattice constant S X1 t denotes discrete time with step size St- We usually set both S x 
and St to unity. The subscript a is an index for velocity; as illustrated in Fig. 1, for the FHP model, a runs 
from 1 to 6. In the figure, collisions that can be obtained by rotations of the input and output states are not 
shown. Also note that momentum must be conserved during collisions, therefore many output states must 
be equal to their input states. These trivial collisions are not shown in the figure. After colliding, particles 
advect to the next site according to their velocities. Fig. 2 illustrates the evolution of the system in one time 
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Input State Output State Input State Output State 



Fig. 1 . Collisions of FHP LG A model. Note that the figure does not included those collisions that can be obtained by 
applying rotations of multiple 7r/3 to input and output states simultaneously. 



Fig. 2. Evolution of FHP LG A model. Solid and hollow arrows represent particles with velocities corresponding to times 
t and t + 1, respectively. That is, the hollow arrows are the final configurations of the initial configurations of solid arrows 
after one cycle of collision and advection. 

step from t to t + 8t. In this figure, solid and hollow arrows represent particles with corresponding velocity 
at time t and t + <^, respectively. The system evolves by iteration of the collision and advection processes. 
According to the collision rules prescribed in Fig. 1, the collision operator, (7 a , can be written as follows: 

(2.2) C a ({n a (x, t)}) = ^2(s' a - s a )£ ss - U<'(1 , 

S,S‘ a 

where s = {si, $ 2 , • • • , $6} and s' = {s^, . . . , s' 6 } are possible incoming and outgoing configurations at 

a given site x and time £, respectively; £ 88 f is a Boolean random number in space and time that determines 
the transition between state s and s f satisfying the following normalization condition: 

(2.3) =1> Vs ‘ 

S' 

The Boolean random number ^ 88 > must also have rotational symmetry, ie., for any states s and s', ^ 88 > is 
invariant if states s and s' are both subjected to simultaneous clockwise or counterclockwise rotations. It is 
obvious that for Boolean numbers n a and s ai the following equation holds: 

(2.4) n s / (1 - n,)* 1 - 8 ’) = b n<jS<J , 

where S n<r8g . is the Kronecker delta symbol with two indices. Therefore, Eq. (2.2) can be written as 
(2*5) ({^a («*m ^)}) = ^ ^ (^a — ^a) £ss f $ns ? 

S,S' 
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Table 2.1 

Collision table for 6-Bit FHP model. 


Input State 

Output State 


010010 

001001 



100100 

010101 

101010 

001011 

100110 


110110 

011011 



101101 


where 8n S = S ni8l S n282 • "S nb8b . Eq. (2.2), or (2.5), is rather abstract. A specific example of the collision 
operator for the two-body collision is 

(2.6) ^ n a +i7l a +47l a n a -\-2'Tta+3'ft'a-\-5 "I" ?i CK -j_2?l'a+5^'a^a+l^'Q:-l-3^'a+4 

(£# “I" £l ) ^a^a+3^a+l^a+2^a+4^a+5 ? 

where n a = 1 — n a is the complement of n a , and £^ 2) are Boolean random numbers which determine 
the outcome of head-on two-body collisions. The Boolean random numbers reflect the randomness of the 
outcomes of the two-body collision. Obviously, for the collision operator to satisfy the complete lattice 
symmetry group statistically (on average), they must satisfy 

(2-7) ® = <C>, 

where (•) denotes the ensemble average. The conservation laws of the particle number, momentum, and 
energy of the LGA micro-dynamics can be written as follows: 

(2.8a) 5>«-Sa) = °, 

a 

(2.8b) ]T(s' a -s a )e a =0, 

a 

(2.8c) ^(s'a -S <*)§ e a =°- 

a 

In practice, the collision can be implemented with various algorithms. One can either use logical op- 
erations [as indicated by Eq. (2.6)], or by table-lookup. The collision rules shown in Fig. 1 can also be 
represented by a collision table, as shown by Table 2.1. In Table 2.1, each bit in a binary number represents 
a particle number n a , a — 1, 2, . . . , 6, from right to left. The limitation of table lookup is the size of the 
table, which is 2 & , where b is the number of bits of the model. Both logic operations and table lookup can 
be extremely fast on digital computers, and especially so on dedicated computers [49, 2]. 

2.2. Hydrodynamics of Lattice Gas Automata. The ensemble average of Eq. (2.1) leads to a 
lattice Boltzmann equation: 

(2.9) fa(x + e a <^, t + St) — fa( x i t) + > 
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where f a = m(n a ) and = {C a ), where m is the particle mass. In addition, correlations among colliding 
particles are assumed to be negligible, «.e., 


( 2 . 10 ) 


(n a np •••%) = (n a )(np) ■ ■ • (n 7 ) . 


The above approximation is equivalent to the celebrated molecular chaos assumption of Boltzmann ( Stosszahlansatz ). 
With the molecular chaos approximation, the lattice Boltzmann collision operator is given by 

(2.11) n«({/«(®, *)}) = 

S, S' a 

where A 88 > = is the transition probability from state s and s'. The hydrodynamic moments are given 
by: 

(2.12) P = £/a, P u - S e ^ a ’ p£ = E ’ 

a a a 

where p, w, and £ are the mass density, the velocity, and the internal energy density, respectively. 

Eq. (2.9) can be expanded in a Taylor series of St up to the second order: 


(2-13) 


(dt + eoC V)/a + ( dt + Ca* V) 2 / a — Cl a 


The equilibrium distribution, fa \ which is the solution of O a ({/ a }) = 0, must be a Fermi-Dirac distribution 
because the system is a binary one [54], that is, 


(2.14) 


fi 0) 


1 

1 + exp(a + b u • e a ) ’ 


where coefficients a and b are functions of p and u 2 in general. Because the coefficients a and b in cannot 
be determined exactly [54], must be expanded in a Taylor series of u — the small velocity (low Mach 
number) expansion. With the small velocity expansion of the equilibrium and through Chapman-Enskog 
analysis, one can obtain the following hydrodynamic equations from the FHP LGA model [19, 54]: 


(2.15a) dtp + V' pu = 0 , 

(2.15b) dtpu + V(gpuu) = —VP + vV 2 pu 1 


where g is a function of p, 
(2.16a) 

(2.16b) 

(2.16c) 


P = c 2 8 p 


u* 

1 - 0-2 

& 


c * = 7i c ’ 

_i /i 


+ 


cS x > 


4 VA ' 2, 

c 8 is the sound speed, c — S x / St, and A is an eigenvalue of the linearized collision operator [38]: 

<9H a I 


Jap — 


dfp 


h=tf 


The defects of the LGA hydrodynamics are obvious from the above equations [Eqs. (2.15a) and (2.16a)]: 
1. Simulations are intrinsically noisy because of the large fluctuation in n a ; 
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2. The factor g(p) is not unity, thus Galilean invariance is destroyed; 

3. Increasing the Reynolds number Re is difficult; 

4. The equation of state depends on u 2 (unphysical); 

5. Spurious (unphysical) conserved quantities exist due to the simple symmetry of the lattice-gas au- 
tomata. 

However, all these defects can be fixed by using more sophisticated LG A models [12], or the other alternative 
— the lattice Boltzmann equation. 


2.3. Lattice Boltzmann Equation. Historically, models of the lattice Boltzmann equation [32, 10] 
evolved from their Boolean counterparts: the lattice-gas automata [19, 54]. Recently, the lattice Boltzmann 
equation has been shown to be a special discretized form of the continuous Boltzmann equation [23, 24]. 

For the sake of simplicity but without lose of generality, the Boltzmann equation with the Bhatnagar- 
Gross-Krook (BGK) approximation [5] is used in the following analysis. The BGK Boltzmann equation can 
be written in the form of an ordinary differential equation: 

(2-17) Z) t / + i/ = I/( ( >), 


where D* = c^+4* V is the Lagrangian derivative along the microscopic velocity / = f(x, t) is the single 
particle distribution function, A is the relaxation time due to collision, and /(°) is the Maxwell-Boltzma.nn 
distribution function: 


(2.18) 


/( °> 


p 

(2ttRT) d / 2 


exp 


20 


in which D is the dimension of the space; 0 = ksT/m is the normalized temperature; T, ks and m are 
temperature, the Boltzmann constant, and particle mass, respectively. The macroscopic variables are the 
moments of the distribution function / with respect to velocity 

( 2 . 19 ) P= J fd£, pu = jtfd£, p$ = ^ J \(-u) 2 fd(. 

Equation (2.17) can be formally integrated over a time interval <^: 

(2.20) f(x + £5 t , 4, t + d t ) = e~ St/x f(x, t) + ^e~ St/x J e t>/x f {0 \x + $t', 4, t + t')dt' . 

Assuming that St is small enough and is smooth enough locally, and neglecting the terms of the order 
0(St) or smaller in the Taylor expansion of the right hand side of Eq. (2.20), we obtain 

(2.21) f(x + (St, 4, t + St)- f{x , 4, t) = --[fix, 4, t) - /<°>(*, 4, t)] , 

T 

where r = A / St is the dimensionless relaxation time. 

The equilibrium distribution function /(°) can be expanded as a Taylor series in w. By retaining the 
Taylor expansion up to tt 2 , we obtain: 

(^tt) oe^) 2 _ 

e 2 < 9 2 20 ■ 

For the purpose of deriving the Navier-Stokes equations, the above second-order expansion is sufficient. 

To derive the Navier-Stokes equations, the following moment integral must be evaluated exactly: 

(2.23) J /0 q) d( , 


(2 ‘ 22) /(eq) i2ir0) D / 2 eXP ( 20/ 


7 



2 


6 


3 


7 


X 

z 

2 

\ 


5 


1 


8 


Fig. 3. Discrete velocities of the 9-velocity model on a square lattice. 


where 0 < m < 3 for isothermal models. The above integral contains the following integral which can be 
evaluated by Gaussian-type quadrature: 

(2.24) I = ( = ex P(-£o/20) , 

^ a 

where 0(4) is a polynomial in and W a and 4 a are the weights and the abscissas (or discrete velocities) of 
the quadrature, respectively. Accordingly, the hydrodynamic moments of Eqs. (2.19) can be computed by 
quadrature as well: 

(2.25) p=J2fa, pu = fa , pO = ~ ~ U f > 

a a ot 

where f a = f a (x 1 1) = W a f(x , £ a , t). We shall use the 9- velocity (9-bit) isothermal LBE model on square 
lattice space as an example to illustrate the derivation of LBE models: the evolution equation (2.21) on a 
discretized phase space and time with a proper equilibrium distribution function leads to the Navier-Stokes 
equations. 

To derive the 9-velocity LBE model, a Cartesian coordinate system is used, and accordingly, we set 
'0(4) = 4™ 4™. The integral of Eq. (2.24) becomes: 

(2.26) I = (V20)( m+n+ V I m I n , 


where 

/*+°° 2 

(2.27) I m = e^TdC, 

J —OO 

and C = 6c/\/20 or ^y/^/20. Naturally, the third-order Hermite formula is the optimal choice to evaluate I m 
for the purpose of deriving the 9- velocity LBE model, ie., I m = Yfj = i • The three abscissas (Q) and 
the corresponding weights (coj) of the quadrature are: 


(228) Cl = -v^72> C2 = 0, C3 = v / 372, 

Wi=v / ir/6, o> 2 = 2 v / 7F/3 , u> 3 = y/n/6. 

Then, the integral of Eq. (2.26) becomes: 

4 8 

(2.29) I = 20 [ui%ip(0) + E E U ^a)] , 


a=l 


a=5 


where 4 a is the zero velocity vector for a = 0, the vectors of a/30 (±1, 0) and a/30 (0, ±1) for a = 1-4, and 
the vectors of a/30(±1, ±1) for a = 5-8. Note that the above quadrature is exact for (m + n) < 5. 

Now the momentum space is discretized with nine discrete velocities {4 a l a = 0, 1, • • • , 8}, as shown in 
Fig. 3. To obtain the 9- velocity model, the configuration space is discretized accordingly, ie., it is discretized 
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into a square lattice with lattice constant S x = a / 3 05t- It should be stressed that the temperature 9 has no 
physical significance here because we are only dealing with an isothermal model. We can therefore choose S x 
to be a fundamental quantity instead, thus V39 = c = S x /St^ or 9 = (? 8 = c 2 ^, where c 8 is the sound speed 
of the model. 

By comparing Eqs. (2.24) and (2.29), we can identify the weights defined in Eq. (2.24): 


(2.30) 

w a 

= 27 T$ exp(£a/2 0)w a , 

where 


f 4/9, 

a = 0, 

(2.31) 

w a = < 

, 1 / 9 , 

a = 1, 2, 3, 4, 



l 1/36, 

a = 5, 6, 7, 8. 


Then, the equilibrium distribution function of the 9- velocity model is: 


(2.32) 

where 


fi eq) = w a f^(x, s a ,t) 


w a p 


' 1+ 3(e a ■ u) + 9(e a ■ uf 


2c 4 


3 u 2 1 
2? J ’ 


(2.33) 


^ a 


' ( 0 , 0 ), 

< (cos sin</> a ) c, 

, (cos(/> a , sin</> a ) a/2c, 


a = 0 , 

a = 1, 2, 3, 4, 
a = 5, 6, 7, 8, 


and (j) a = (a — l)7r/2 for a = 1-4, and (a — 5)7r/2 + 7r/4 for a = 5-8, as shown in Fig. 3. The Navier-Stokes 
equation derived from the above LBE model is: 


(2.34) pdtu + pu-Vu = — Vp + pvV 2 u , 

where the equation of state is the ideal gas one, p = c^p, the sound speed c 8 — c/a/ 3 ? and the viscosity 
(2-35) v = | (t - 1/2 )c 8 x = {t- 1 / 2 )c; % 

for the 9-velocity model. 

Similarly, we can also derive two-dimensional 6-velocity, 7-velocity, and three-dimensional 27-velocity 
LBE models [24]. 

In the above derivation, the discretization of phase space is accomplished by discretizing the momentum 
space in such a way that a lattice structure in configuration space is simultaneously obtained. That is, the 
discretization of configuration space is determined by that of momentum space. Of course, the discretiza- 
tion of momentum space and configuration space can be done independently. This consideration has two 
immediate consequences: arbitrary mesh grids and significant enhancement of the Reynolds number in LBE 
hydrodynamic simulations. 

To implement arbitrary mesh grids with the LBE method, one first discretizes the configuration space 
by generating a mesh adapted to the physics of the particular problem. Then at each grid point, one can 
discretize momentum space as before. Now, a local LBE is built on each mesh grid point. The evolution of this 
discretized Boltzmann equation (DBE) consists of the following three steps. The first two steps are the usual 
collision and advection process as in the previous LBE models. After collision and advection, interpolation 
follows. The interpolation process is what distinguishes the DBE from the LBE method. Because the mesh 
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grids can be arbitrary, the distribution function f a at one mesh grid point, say a?*, cannot go to another grid 
point in general through the advection process as it can in previous LBE models. Therefore, the interpolation 
step becomes necessary to construct t) on each and every mesh grid point from f a (xi + e a St , t) after 

the advection process. Of course, interpolation brings in additional numerical error, but it can be justified so 
long as the error induced by interpolation does not affect the DBE algorithm as a whole [25]. In addition, the 
separate discretization of momentum and configuration space allows one to increase the Reynolds number 
significantly in numerical simulations without enlarging mesh sizes or decreasing the viscosity by adjusting 
r [25]. In other words, the limitation posed by the lattice Reynolds number is completely overcome [26] and 
the stability of the LBE method is greatly improved [25, 26]. 

2.4. Boundary Conditions and Turbulence Modeling. The boundary conditions for the lattice 
gas automata are relatively simple. The simplest boundary conditions for the velocity field are bounce 
back and reflective boundary conditions. The former mimics the no-slip boundary condition and the latter 
simulates the slip boundary condition. In the bounce back boundary condition, a particle hitting a wall 
simply reverses its momentum, whereas in the reflective boundary condition the particle hitting a solid 
boundary only reverses its momentum normal to the boundary and maintains its momentum tangential to 
the boundary. In addition, fully developed boundary conditions for the velocity field can be obtained by 
equating the distribution functions in two adjacent lattice lines. The pressure boundary condition can be 
implemented via the ideal-gas equation of state by setting appropriate mass densities in the desired locations. 
Body forces can also be implemented by probabilistically making particles gain momentum along the forcing 
direction. The Boolean nature of the lattice gas automata sometimes makes the boundary conditions difficult 
to alter. Furthermore, the accuracy of the aforementioned boundary conditions could be of first order (in 
space), instead of second order. 

The lattice Boltzmann method has greater freedom than lattice gas automata methods to implement 
various boundary conditions accurately [33]. Various conventional techniques such as interpolation or ex- 
trapolation can be applied to improve the accuracy of the LBE boundary conditions [33]. Furthermore, 
arbitrary nonuniform meshes [25] and mesh refinement [18] can be implemented in the LBE method. 

There are two types of turbulence models implemented in the lattice Boltzmann method: subgrid mod- 
eling [27] and wall function modeling [3]. The subgrid modeling is illustrated as follows by using the example 
of the Smagorinsky subgrid model. In a subgrid model, the total viscosity v can be divided in two parts, 

(2.36) v - v 0 + v s , 

where uq and v 8 are the physical viscosity (for the resolved scales) and eddy viscosity (for the unresolved 
subscales), respectively. The eddy viscosity is obtained as follows: 

(2.37a) v g = C 2 5l\S\, 

(2.37b) |5| = 

(2.37c) §ij = — ( djiii + ditij) , 

where u is the resolved velocity field, §ij is the corresponding strain-rate tensor, and C > 0 is the Smagorinsky 
constant. The relaxation time is evaluated at each grid as follows: 

(2.38) r = 6(u 0 + v 8 ) + - , 

where uq is given a priori , as determined by the Reynolds number of the simulated flow. The relaxation 
time now varies locally to model the subgrid-scale dynamics, depending on the flow field. 
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An alternative to the above subgrid modeling of turbulence is through a wall function [3]. Instead of 
resolving the boundary layer adjacent to the wall, the computation starts at the first cell above the surface, 
where a “slip” velocity U 8 is typically nonzero [45]. The turbulent law of the wall is assumed to hold at the 
location y 8 of the first cell off the wall; hence, the shear stress at y 8 is the same as the wall shear stress r w . 
A local skin-friction coefficient C'j is defined through the relation 

(2.39) t w = cf-pUl 

The local friction velocity U r is defined as 


(2.40) U T = ^ 

Application of the law of the wall at the distance y 8 from the wall yields the logarithmic profile of the 
boundary layer 


(2.41) 


U T 


I In (TlKi ] + Bk I 


where k is the von Karman constant (k ss 0.4) and B is a constant of integration (B « 5.0). From Eqs. 
(2.39) and (2.40), 


(2.42) 


U, 

U T 


and the friction coefficient can be expressed as 



(2.43) 


Cf 


2 K? 

In(^) +Bk 


In a simple mixing-length model a turbulent, or eddy viscosity Ut is proportional to the mean velocity gradient 
multiplied by the square of a length scale. If y 8 is in the log layer, ny 8 is an appropriate length scale and 
U r /(Ky 8 ) is the mean velocity gradient, hence 


(2.44) 


Ut — Ky 8 U T . 


Anagnost et al. [3] assume that the lattice viscosity u* is equal to the eddy viscosity u t at the height y 8 of 
the first grid cell. Why such an assumption is necessarily valid is not clear; however, this assumption leads 
to the basic form of the relationship for the local skin-friction coefficient 


(2.45) 


C'f 


2k? 

m£)+£«t 


The effect of local adverse pressure gradients can be incorporated into the model by modifying the relationship 
between the lattice viscosity and the eddy viscosity. The modified equation for C'j is 


(2.46) 


C'} 


2 k? 

P n {(^) (! + x(^p))} +Bk} 2 ' 


where the function x{®xP) is negative in an adverse pressure gradient and zero otherwise. The C’j effectively 
determines a “slip” at the boundary. This “slip” boundary condition can be accomplished by a proper ratio 
between bounce-back and reflection of the distribution function colliding with the wall. 
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3. PowerFLOW. PowerFLOW is a software package based on lattice gas cellular automata [17]. Pow- 
erFLOW uses 3 speeds in three-dimensional space [45]. The equilibrium distribution function in the software 
is of Maxwellian type: 

(3-1) fi e9) = pvJa(T) exp (T( 2e a ■ u - u 2 )^ , 


where w a (T ) can be determined by the conservation constraints [12]. In particular, for a 3 speed model [45]: 


(3.2) 


w a (T) 


' £(3T 2 -3T+1), 

< ^2T(2-3T), 

. &r(3T-l), 


|e a | = 0, 

|e a | = 1, 
\e a \ = V% 


where do Is an adjustable positive constant and is usually set to 6. The distribution function in PowerFLOW 
is represented by an eight-bit integer [45]. Some special care must be taken to ensure exact preservation of 
the conservation laws. 

PowerFLOW uses a structured, cubic Cartesian mesh (voxels), and a grid refinement scheme which 
refines the grid size in each spatial dimension equally and successively by a factor of 2. The boundary 
surface is tessellated into triangular planes (surfels), of which each side has a length that is less than 1/10 
the linear dimension of a voxel in the region. 

Various boundary conditions are implemented in PowerFLOW by specifying a velocity profile, static 
pressure, and mass flux at the boundary [17]. At the boundary, momentum fluxes are computed and 
redistributed to enforce desirable boundary conditions [13]. Because PowerFLOW is a commercial software 
package, and source code was unavailable to the authors, specific implementation details discussed in Sections 
2 and 3 may not be completely relevant to the release of the PowerFLOW software used in this work. 


4. CFL3D. A standard computer code, CFL3D [48, 51, 47, 40, 39, 29], was used to compute flow 
fields for comparison with PowerFLOW results. The computer code CFL3D was developed at NASA Lan- 
gley Research Center to solve the three-dimensional, time-dependent, thin-layer (in each coordinate direc- 
tion) Reynolds-averaged Navier-Stokes (RANS) equations using a finite- volume formulation. The code uses 
upwind-biased spatial differencing for the inviscid terms and flux-limiting to obtain smooth solutions in the 
vicinity of shock waves. The viscous derivatives are computed by second-order central differencing. Fluxes 
at the cell faces are calculated by Roe’s flux-difference-splitting method [36]. An implicit three-factor ap- 
proximate factorization method is used to advance the solution in time. Here, only steady solutions have 
been computed and convergence is accelerated with local time-stepping, grid sequencing, and multigridding. 


4.1. Turbulence Modeling. The CFL3D flow solver provides a wide range of turbulence models from 
which to choose. In this work, most of the cases employed the model of Spalart and Allmaras [42]. One case 
used Wilcox’s k-cv model [52]. 

In the model of Spalart and Allmaras [42] the turbulent eddy viscosity is obtained directly from the 
solution of a transport equation. The model is tuned for aerodynamic flows, is robust, and is applicable 
to flows with complex geometry. Several studies for both steady and unsteady flows (Rogers et al. [37] 
and Rumsey et al [40]) show the model to have better predictive capability than algebraic and other one- 
equation models for a wide range of aerodynamic flows. The CFL3D implementation uses Version la of the 
model [43]. 

The Wilcox k-oo model [52] is a two-equation model, one for the turbulent kinetic energy fc, the other 
for which is proportional to the turbulent energy dissipation rate per unit kinetic energy in the CFL3D 
implementation [29]. The resultant turbulent eddy viscosity is the ratio k/oj. 
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Although not required by CFL3D, in cases in which a turbulence model is used, the flow is considered 
fully turbulent from the onset, ie., no laminar region is assumed. In addition, both the Spalart-Allmaras 
model and the Wilcox model are integrated all the way to the wall. As with all turbulence models in CFL3D, 
the solution of either turbulence model is decoupled from the solution of the flow equations. Information 
from the flow equations that is required in the turbulence equation lags one iteration. This is not a problem 
for sufficiently converged solutions. Additional implementation details are available in the CFL3D manual 
[29]. 

4.2. Boundary Conditions. As shown in Fig. 5, a C-grid was used for the CFL3D calculations. This 
grid requires boundary conditions to be defined on three distinct geometrical surfaces: (1) the wing body, (2) 
the downstream plane, and (3) the external surface that defines the upstream, and upper and lower portions 
of the computational domain. 

For the inviscid calculations with CFL3D, slip wall conditions (CFL3D boundary-condition type 1005) 
are used on the wing body (surface 1). These conditions require that the component of velocity normal to 
the wall be equal to zero. In CFL3D this boundary condition approximates the cell-face boundary value for 
the density as the density of the nearest cell-center point on the grid. The boundary values for the pressure 
are obtained similarly. For the inviscid cases, the external computational surfaces (surfaces 2 and 3) are 
modeled with the use of local one-dimensional Riemann invariants (CFL3D boundary condition type 1003). 
The conditions outside the computational domain are taken from the prescribed free-stream Mach number 
and flow direction. 

For the viscous calculations, no-slip, no-penetration, adiabatic wall conditions (CFL3D boundary-condition 
type 2004 with twtype = 0 and cq — 0) were used on the wing body (surface 1). In the viscous cases, the 
condition imposed at the downstream plane (surface 2) involved a zeroth-order extrapolation from the com- 
putational interior (CFL3D boundary-condition type 1002). The far-field, upstream computational surface 
(surface 3) had its boundary condition modeled with the local one-dimensional Riemann invariant, as with 
the inviscid cases. 

More details regarding the implementation of the boundary conditions can be found in the CFL3D 
manual [29]. 

5. Two-Dimensional Simulation of Flow Past NACA-0012 Airfoil. 

5.1. NACA-0012 Airfoil. The NACA-0012 airfoil is a popular wing section that is used for many 
purposes. Coordinates and a generation formula are given in Ref. [1]. Figure 4 shows a NACA-0012 wing 
section. The NACA-0012 airfoil has zero camber and a maximum thickness to chord ratio of 12 percent. A 12 
percent thickness ratio results in maximum lift for many 4 and 5 digit NACA wing sections [1]. Experimental 
lift versus angle-of-attack curves, reproduced in graphical form in Ref. [1], suggest that on aerodynamically 
smooth wings, the NACA-0012 airfoil does not experience massively separated flow up to an angle of attack 
of 16° for Reynolds numbers of 3, 6, and 9 million. For the wing with a standard roughness, massive flow 
separation occurs at an angle of attack of about 12° at a Reynolds number of 6 million. Because the NACA- 
0012 airfoil is symmetric, it generates no lift at zero angle of attack. At an angle of attack of 3°, the lift 
coefficient is expected to be about 0.35. At an angle of attack of 7°, the lift coefficient is expected to be 
between about 0.74 and 0.78, depending upon conditions. 

5.2. Computational Domain and Boundary Conditions. A series of 2-D simulations of flow past 
the NACA-0012 airfoil at low speed were performed by using PowerFLOW and the results are compared 
with those from CFL3D. The Mach number M is set to 0.2 in the simulations. Several values of angle of 
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Fig. 4. NACA-0012 airfoil. 


attack, a, are selected in the simulations: 0°, 3°, 7°, and 12°. 

Both CFL3D and PowerFLOW are three-dimensional codes. However, the simulations of the flow past 
the airfoil are two-dimensional. In the case of CFL3D, this is simply done by using the two-dimensional 
solver option in the code. In the case of PowerFLOW, the two-dimensional simulations were accomplished by 
shrinking the third dimension to two cells and employing a periodic boundary condition [46]. The simulations 
were done on an SGI Origin2000 with 16 R10000 processors. The CPU times and resolutions used in the 
simulations of both CFL3D and PowerFLOW are summarized in Table 5.1. 

For PowerFLOW, a typical computational domain is L x x L y = 9 x 8 in units of chord length. Pow- 
erFLOW uses Cartesian grids. The resolution at the surface is approximately 500 points per chord. In the 
PowerFLOW simulations there are seven variable resolution regions, i.e., the grid spacing is varied from 2° 
to 2 -6 . The resolution used in PowerFLOW is usually 383,724 voxels (grid points) and 1,275 surfels (surface 
elements). The aspect ratio of grid spacings is always 1 : 1 in the PowerFLOW simulations. 

The boundary conditions used in PowerFLOW are: a free-slip boundary condition for upper and lower 
boundaries; a driving boundary condition with a fixed velocity at the inlet boundary, and a fixed static 
pressure for exit boundary condition. At the low Reynolds number (Re = 500), a no-slip boundary condition 
is used for the airfoil surface. At high Reynolds numbers (Re > 0.5 x 10 6 ), the “turbulent wall boundary 
condition” using a wall function is applied to the airfoil surface. 

PowerFLOW is a time explicit scheme with a small time step: the Courant-Friedrichs-Lewy number is 
one with respect to the grid. PowerFLOW has significant fluctuations in its instantaneous results because 
of its integer nature. Therefore hydrodynamic outputs of PowerFLOW are usually averaged over a number 
of time steps. The interval for the time-average for the NACA-0012 airfoil simulations is 1,000 steps unless 
otherwise specified. 

Figure 5 shows the body-fitted computational mesh (C-grid) used in CFL3D simulations. With the grid 
nondimensionalized with the chord length, the computational domain consists of a half disc of radius 20 in 
front of the airfoil and a rectangle of 20 x 40, as shown in Fig. 5. The resolution is typically 373 x 141, 
or 52,593 cells. The grid spacing normal, t/i, and tangential, aq, to the surface near the middle of chord 
are approximately 0.005 and 0.015 for in viscid cases, respectively, and 0.00012 and 0.008 for viscous cases, 
respectively. With these fine grids, CFL3D can resolve the boundary layer in the simulations. 

Although all flows possess some unsteadiness, for attached flows the fluctuations typically are on a time 
scale that is so small that they are usually of no practical interest for aerodynamic quantities like lift and 
drag. Most experimental force data are actually time-averaged quantities. Some sensors implicitly average 
because they are incapable of resolving the small, unsteady fluctuations. To obtain mean data from a 
time-accurate CFD simulation requires performing a similar time-average on the numerical results. This is 
the approach used in PowerFLOW. Although only fractions of a second are needed to obtain a significant 
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Fig. 5. A typical body -fitted mesh used in CFL3D. (left) the entire computational domain with grid resolution 373 x 14 1, 
and (right) the mesh around the airfoil , which is the center part of the domain. 


sample in an experiment, considerable computer resources are required to generate all of the required samples 
from a numerical simulation. To simplify and accelerate the solution procedure, CFD schemes commonly 
solve equations that have been time-averaged so that all of the dependent variables are actually time-mean 
quantities. In the case of the Navier-Stokes equations applied to turbulent flows, this process is called 
Reynolds averaging and produces several additional terms which must be modeled to close the system of 
equations. This modeling of the effects of the small scale, unsteady turbulent flow on the mean quantities 
is nontrivial, but many models have been developed and validated for many important flows. Because the 
equations are nonlinear, the typical solution procedure is to iterate until the solution becomes steady. Many 
acceleration techniques have been developed that take advantage of the fact that the intermediate solution 
need not be time accurate. This is the approach used in CFL3D. 

5.3. Simulation Results. The local results discussed below include u and t;, the x and y components 
of the velocity field normalized by the inlet velocity Uoo\ and the pressure coefficient 


(5.1) 




p- Poo 

-O U 2 ’ 

2 Voodoo 


where p is the local pressure, p^ is the far-field pressure, and p^ is the far-field density. The velocity and 
pressure are measured at five stations: x/L — 0, 1/4, 1/2, 3/4, and 1, where L is the chord length. In 
addition, the lift coefficient 


(5.2) 


and the drag coefficient 
(5.3) 


Fit A 


Fd/A 

hPooU*,’ 


are obtained by integrating the lift force F} and the drag force per unit span over the surface of the airfoil. 

The resultant lift coefficients Ci and drag coefficients Cd are summarized in Table 5.1. The Reynolds 
numbers Re in these simulations range from 500 to 0.5 x 10 9 . For Re = 500 (Case la and lb in Table 5.1), 
the results were obtained by directly simulating the laminar flow. For Re > 500 x 10 3 , turbulence modeling 
was used in both PowerFLOW and CFL3D (Cases 2 - 6 in Table 5.1). PowerFLOW uses wall functions, 
described in Sec. 2.4, to model the boundary layer. CFL3D incorporates various turbulence models in the 
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Table 5.1 

Summary of results for simulations for flow past a NACA-0012 airfoil. (Top) Drag coefficient C d and lift coefficient Ci 
obtained from CFL3D and PowerFLOW. The numbers in parentheses are bounds estimated from the unsteady results of the 
PowerFLOW simulations. (Bottom) Resolutions and CPU times of simulations using CFL3D and PowerFLOW. 


Case 

a 

Re 

CFL3D 

PowerFLOW 

c d 

Cl 

c d 

Ci 

la 

0° 

500 

0.17618 

0.115X10" 6 

0.171721 

0.22703xl0- 3 

lb 

0° 

500 

0.1741 

-0.538X10" 5 

0.180716 

— 0.211xl0 -3 

2a 

0° 

0.5 x 10 6 

0.016268 

0.1769 x 10-® 

0.0212145 

-0.2 x lO -3 

2b 

0° 

0.5 x 10 6 

0.012174 

0.9241 x 10" 9 

0.0201036 

-0.3518 x lO -3 

2c 

0° 

0.5 x 10 6 

0.012815 

-0.616 x 10" 7 

— 

— 

3 

3° 

0.5 x 10 6 

0.012986 

0.32372 

0.0132659 

0.188275 

4a 

3° 

in viscid 

0.001256 

0.3608 

— 

— 

4b 

3° 

0.5 x 10 9 

— 

— 

-0.014939 

0.351683 

4c 

3° 

0.5 x 10 9 

— 

— 

0.00304585 

0.0512701 

5 

7° 

0.5 x 10 6 

0.015739 

0.7449 

(-0.017, 0.0019) 

(0.72, 1.2) 

6 

12° 

0.5 x 10 6 

0.027468 

1.1809 

(-0.054, -0.12) 

(0.7, 1.56) 

Case 

a 

Re 

CFL3D 

PowerFLOW 

N e x N r 

CPU Time s (hrs) 

Voxel xSurfel 

CPU Time s (hrs) 

la 

0° 

500 

257x65 

3500 (1.0) 

159060x828 

36795 (10.2) 

lb 

0° 

500 

373x141 

23143 (6.4) 

418800x1275 

70900 (19.7) 

2a 

0° 

0.5 x 10 6 

257x65 

1650 (0.5) 

145748x828 

23003 (6.4) 

2b 

0° 

0.5 x 10 6 

373x141 

17671 (4.9) 

383724x1275 

63034 (17.5) 

2c 

0° 

0.5 x 10 6 

373x141 

17662 (4.9) 

— 

— 

3 

3° 

0.5 x 10 6 

373x141 

29962 (8.3) 

383724x1275 

60772 (16.9) 

4a 

3° 

inviscid 

225x65 

1630 (0.5) 

— 

— 

4b 

3° 

0.5 x 10 9 

— 

— 

383724x1275 

66940 (18.6) 

4c 

3° 

0.5 x 10 9 

— 

— 

63612x604 

112528 (31.3) 

5 

7° 

0.5 x 10 6 

373x141 

29962 (8.3) 

383724x1275 

61627 (17.1) 

6 

12° 

0.5 x 10 6 

373x141 

22680 (6.3) 

383724x1275 

67320 (18.7) 


code. Two turbulence models were used in the simulations of flow past a NACA-0012 airfoil: the Spalart- 
Allmaras model and the k-oo model. The k-oo model was only used in Case 2c of Table 5.1 to demonstrate the 
level of variability in the CFL3D result with different turbulence models. The results show little dependence 
on the turbulence model for this flow. The remainder of CFL3D simulations used the Spalart-Allmaras 
model. 

Case 1: Re = 500, a — 0. The CFL3D simulations used two resolutions: 257 x 65 and 373 x 141, and 
PowerFLOW used 159,060 voxels and 828 surfels (Case la), and 418,800 voxels and 1275 surfels (Case lb). 
The value of drag coefficient Cd obtained from PowerFLOW is about 3.8% larger than the value obtained 
from CFL3D. Although a symmetric airfoil at zero angle of attack should generate zero lift, the PowerFLOW 
lift coefficient varied between ±0.0002. The CFL3D solutions maintained zero lift to an accuracy almost two 
orders of magnitude better. 
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Fig. 6. Time history of Cd and Ci in the PowerFLOW simulation. Re = 500, M = 0.2, and a = 0°. 

At a Reynolds number Re = 500, the flow is laminar and the results were obtained by direct numerical 
simulation — no turbulence modeling was used. Along the surface of airfoil, the no-slip boundary condition 
is used in PowerFLOW. After 25,000 time steps, a steady-state solution was attained. Figure 6 shows the 
time history of Cd and C*. The variations of Cd and Ci are on the order of 10 -4 . 

Figures 7, 8, and 9 show u(x , $/), v(x , t/), and C p (x 1 y) at five vertical sections of x/L — 0, 1/4, 1/2, 
3/4, and 1. These figures show the results with the higher resolution case (Case lb). The differences 
between the results of PowerFLOW and CFL3D are already visible at Re = 500. The difference in the x- 
component of velocity u is relatively small. However, the relative difference in the y-component of velocity v 
is significantly greater, and the difference is amplified downstream. Similar discrepancies exist in the vertical 
pressure distribution and the surface pressure distribution. Although the no-slip condition is satisfied in the 
CFL3D simulations, the figures do not show u — > 0 near the airfoil surface. This is because the data were 
extracted at uniformly spaced points in the y direction, and a point on the airfoil surface is usually not 
included. This is also true for the PowerFLOW data. 

Case 2: Re = 0.5 x 10 6 , a = 0°. When Re = 0.5 x 10 6 , turbulence modeling was used. CFL3D used 
the Spalart-Allmaras model (Case 2a and 2b in Table 5.1) and k-uo model (Case 2c). PowerFLOW used the 
wall function discussed in Section 2.4 (Case 2a and 2b). The value of drag Cd obtained from PowerFLOW is 
about 66% larger than the value obtained from CFL3D. The value of Ci from the PowerFLOW simulation 
is 2 to 6 orders of magnitude larger than that of CFL3D. As before, for a = 0°, the value of Ci should be 0. 

Figures 10, 11, and 12 show u(x , t/), v(x , ?/), and C p (x , y) at the five vertical sections, for the results with 
higher resolution (Case 2b). The differences between the CFL3D results and PowerFLOW results are quite 
significant on the surface of the airfoil, and relatively negligible in the flow region away from the surface. 
Obviously, the differences in the flow fields near the boundary result in the significant difference in the value 
of drag Cd * This raises the question of adequacy of turbulence modeling by a wall function in PowerFLOW, 
even for cases in which the flow is attached. 

Case 3: Re = 500 x 10 3 , a — 3°. The results of the drag coefficient Cd obtained from the two methods 
are very close in this case: Cd from PowerFLOW is only about 2.1% larger than that from CFL3D. With 
a = 3°, the lift coefficient is non-zero. The lift coefficient Ci obtained from PowerFLOW simulations is 
about 58% smaller than the result from CFL3D. Because drag includes both a frictional component and a 
component related to lift, the close agreement of the drag coefficient is likely to be fortuitous. 

Figures 13, 14, and 15 show u(x , t/), v(x , ?/), and C p (x , y) at the five vertical sections. The velocity fields 
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obtained from PowerFLOW and CFL3D simulations significantly disagree near the airfoil surface, which in 
turn leads to the disagreement in the region away from the boundary. Furthermore, qualitative disagreement 
in pressure C p (x, y) starts to appear (e.<?., C p (x, y) at x = 3/4 and 1 in Figs. 15). Figure 16 shows the 
streamlines of PowerFLOW and CFL3D simulations in the vicinity of the airfoil. The streamlines do not 
stay completely attached to the airfoil surface in the PowerFLOW simulations. Note that the PowerFLOW 
calculation shows separated flow near the trailing edge while the CFL3D solution remains completely at- 
tached. The separation shown in Fig. 16 could also be surmised from the surface pressure distribution near 
the trailing edge in Figure 15. 

Case 4 : Inviscid flow , Re = oo, a = 3°. In the CFL3D simulation of inviscid flow, an Euler scheme is 
used by neglecting the viscous flux terms. In the PowerFLOW simulation, the Reynolds number is 0.5 x 10 9 , 
and the free slip boundary condition is applied at the surface of the airfoil [17] to mimic the inviscid effect. 
In this case, the value of Cd from one PowerFLOW simulation is negative, as indicated in Table 5.1 (Case 
4b). 

Figures 17, 18, and 19 show u(x , t/), v(x , y ), and C p {x, y) at the five vertical sections, for Case 4a 
(CFL3D) and 4b (PowerFLOW). The results of u(x , y) from CFL3D and the PowerFLOW simulations are 
qualitatively different near the airfoil surface. The CFL3D results were compared against those from other 
CFD codes, and the agreement was excellent. 

Figure 20 shows the time history of Cd and Ci in PowerFLOW simulations with two different resolutions 
and computational domain sizes. In the first case (Case 4b), the resolution is 383,724 voxels and 1,275 
surfels, the resolution on the airfoil surface is 500 points along the chord, the computational domain is 9 x 8, 
and there are seven variable resolution regions (2° to 2~ 6 ). The drag coefficient Cd converges to a negative 
value after about 20,000 iterations, as shown in the left part of Fig. 20. 

In the second case (Case 4c), the resolution is 63,612 voxels and 604 surfels, the resolution on the airfoil 
surface is 200 points along the chord, the computational domain is 21 x 20, and there are seven variable 
resolution regions (2° to 2 -7 ). This case was run to investigate whether the original computational domain 
was too small. With a coarser resolution but a larger computational domain, the results of Cd and Ci do 
not seem to converge after 60,000 iterations, as shown in the right part of Fig. 20. The values of Cd and Ci 
given in Table 5.1 (Case 4c) are the final values of the simulations. They are apparently far away from the 
correct values. It is also noted that the CPU time for the PowerFLOW simulation is almost doubled in this 
case, despite the fact that the resolution is much coarser. 

Case 5: Re = 500 x 10 3 , a = 7°. With a = 7°, the mean flow should be steady and attached. How- 
ever, the PowerFLOW simulation displays massive flow separation, and flow field becomes highly unsteady. 
Figure 21 shows a snapshot of streamlines in both the PowerFLOW and the CFL3D simulation. The two 
results are completely different from each other. Because PowerFLOW failed to converge to a steady mean 
flow after 30,000 time steps, the values of C\ and Cd can only be given in their estimated bounds, as shown 
in Fig. 22. 

Figures 23, 24, and 25 show «(#, t/), v(x, y), and C p (x, y) at the five vertical sections. In this case the 
flow fields obtained from the two methods do not have any resemblance; while the solution from CFL3D 
remains completely attached, the solution from PowerFLOW undergoes a massive separation, as also shown 
in Fig. 21. 

Case 6: Re = 500 x 10 3 , a = 12°. At a = 12°, the observation is similar to the previous case of a = 7°; 
in the PowerFLOW simulation the flow separates massively while it should not. The mean flow fails to 
converge to a steady state, as it should. Therefore, the flow fields obtained in PowerFLOW simulation are 
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qualitatively incorrect, as shown in Figures 26, 27, and 28. 

6. Discussion and Conclusion. In this work we present the results of a thorough investigation using 
PowerFLOW to simulate two-dimensional flow past a NACA-0012 airfoil with a Reynolds number ranging 
from 500 to 0.5 x 10 9 . The PowerFLOW results were compared with the results obtained from CFL3D. 

In all cases we observed deficiencies with the PowerFLOW calculations. We hypothesize that the princi- 
pal culprit for the observed deficiencies is the the lack of wall-normal grid resolution near the surface of the 
airfoil. Because PowerFLOW requires grid cells to be square, properly sizing the grid cells in the wall-normal 
direction requires grid cells in the streamwise direction to be so small as to make even a two-dimensional 
calculation infeasible. 

The lack of wall-normal grid resolution is probably the cause for the near- wall velocity and pressure 
errors that were shown in Figs. 7, 8, and 9. Errors in the near- wall region degrade the accuracy of the overall 
aerodynamic quantities, such as the lift and drag coefficients. Although the use of wall functions for the 
high-Reynolds number turbulent cases relieves some of the resolution requirements for those cases, the data 
suggest that the resolution used was still insufficient. 

Figure 29 shows a velocity profile normal to the wing surface as computed by CFL3D. The variables 
plotted are in wall coordinates (i.e., y + = u T y/v and u + = U/u T where u r — \J r wa ii / p and r wa ii is the wall 
shear stress), which emphasize the profile in the near- wall region. The vertical lines indicate the locations of 
the first through fifth grid points from the surface in the PowerFLOW grid. Although the first PowerFLOW 
grid point is located at an appropriate position for implementing a wall model, the PowerFLOW grid only 
has about 5 points in the entire boundary layer. That level of grid resolution is inadequate for properly 
computing the boundary layer. 

In this flow, incorrectly computing the boundary layer produces catastrophic problems. Flow separation 
occurs when the downstream momentum of the fluid in the boundary layer is insufficient to maintain the 
downstream movement of the fluid against an adverse pressure gradient. The momentum of the fluid in 
the boundary layer is decreased by the shear stress at the wall, but is increased by the (usually turbulent) 
diffusion of momentum from the free stream into the boundary layer. Incorrectly computing the boundary 
layer can dramatically alter where the flow first separates and the degree of separation. After massive flow 
separation occurs, the computed flow has little resemblance to its attached-flow counterpart. 

In the test problems studied here, minor separation was observed in the PowerFLOW results for a case 
with a three-degree angle of attack. Massive separation was observed in the PowerFLOW results for a seven- 
degree angle of attack case. After separation occurred, even the gross aerodynamic quantities could not be 
computed well. 

Two additional tests were carried out at NASA Langley Research Center: (1) simulation of three- 
dimensional cavity flow [28], and (2) three-dimensional flow past the ROBIN configuration (ROtor Body 
INteraction), a generic helicopter configuration consisting of a body and an engine-transmission pylon [34]. 

The geometry of the cavity flow is very simple. The Mach number and the unit Reynolds number of 
the problem are 0.3 and 6 x 10 6 /ft., respectively. The resolution required to obtain reasonable results is 
9,240,000 voxels and 286,500 surfels. A run of 20,000 time steps required approximately 60 hours on an SGI 
Origin2000 with 16 250MHz processors. Although quantitative comparisons with experimental results of 
mean wall pressure distributions, the sound pressure level, and length scales of the flow were not satisfactory 
[28], the simulation did capture the acoustic interaction with the shear layer in the flow. Constraints on time 
and memory size precluded further simulations with finer resolutions. 

The ROBIN is a generic three-dimensional helicopter configuration. The Reynolds number of the flow is 
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9.2 x 10 6 based on the body length (2 meters). Five levels of grid resolution were used in the PowerFLOW 
simulations (2° to 2 -4 ). The resulting grid had 4.3 x 10 6 voxels. The PowerFLOW simulations were compared 
with results obtained by using Navier-Stokes solvers on a structured grid (CFL3D) and an unstructured grid 
(USM3D), and with experimental observations. The main concerns with the PowerFLOW simulations of 
the flow past ROBIN are: (1) too large a value for the maximum pressure coefficient C v \ (2) an excessive 
amount of oscillation in the pressure coefficient distributions and skin friction coefficient distributions; (3) 
significant differences in surface flow patterns and nondimensional stagnation pressure distributions from 
those predicted by the Navier-Stokes solvers; and (4) drag and lift coefficients that differ significantly from 
those predicted by the Navier-Stokes solvers [34]. 

Based on our investigation of the two-dimensional flow past the NACA-0012 airfoil, it is our conclusion 
that PowerFLOW at its present state is inadequate for aerodynamic applications in two aspects: performance 
and accuracy. PowerFLOW is generally much less accurate than CFL3D and much slower as well. 

Some successful examples using PowerFLOW include simulating the flow past a vehicle [3]. The reason 
for this success, we suspect, is that in the case of the flow past a vehicle, the bluff body has geometric 
discontinuities that dictate the location of flow separation. Therefore, the flow separation, which is the most 
important flow feature in this case, is not affected by other factors as much. In this case, PowerFLOW seems 
to be able to deliver better results [3]. However, when dealing with a smooth or streamlined body such as 
an airfoil, PowerFLOW does not provide reliable results. We also observed the fact that simulations using 
PowerFLOW to obtain a steady mean state generally require much larger memory and longer CPU time 
than the CFL3D calculations. 

Although the results of the aerodynamic simulations revealed a less than adequate performance of 
PowerFLOW compared with CFL3D, further advancements in LG A and LBE methodology may enable 
this approach to become a viable alternative to traditional CFD. LBE methods have already been used 
successfully for applications involving porous media and bluff bodies, and further research should enable 
successful applications to streamlined aerodynamic bodies. The current results are presented to illustrate 
the current state of these methods, not to predict their future. Some researchers have already obtained good 
solutions for airfoil flows with Reynolds number Re between 500 and 5000 by using interpolation schemes 
at boundary and multi-block grid to handle surface curvature [50]. Continued refinements by the LBE 
community and Exa corporation will undoubtedly have occurred by the time this paper is disseminated. 
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Fig. 7. NACA-0012 airfoil with angle of attack a = 0°, Re = 500, and M = 0.2. The x -component of velocity u(x y y) 
at x/L = 0.0, 0.25, 0.5, 0.75, and 1.0, where L is the chord length. The resolution is 418,800 voxels and 1,275 surfels for 
PowerFLOW and 373 x 141 for CFL3D. 
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Fig. 8. NACA-0012 airfoil with angle of attack a = 0° , Re = 500, and M = 0.2. The y -component of velocity v(x> y) 
at x/L = 0.0, 0.25, 0.5, 0.75, and 1.0, where L is the chord length. The resolution is 418,800 voxels and 1,275 surfels for 
PowerFLOW and 373 x 141 for CFL3D. 
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Fig. 9. NACA-0012 airfoil with angle of attack a = 0° , Re = 500, and M = 0.2. The pressure coefficient C p (x , y ) at 
x/L = 0.0, 0.25, 0.5, 0.75, and 1.0, where L is the chord length , and along the surface. The resolution is 418,800 voxels and 
1,275 surfels for PowerFLOW and 373 x 141 for CFL3D. 
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Fig. 11. NACA-0012 airfoil with angle of attack a = 0°, Re = 0.5 x 10 6 , and M = 0.2. The y-component of velocity 
v(x, y) at x/L = 0.0, 0.25, 0.5, 0.75, and 1.0, where L is the chord length. The resolution is 383,724 voxels and 1,275 surfels 
for PowerFLOW and 373 x 141 for CFL3D. 










Fig. 12. NACA-0012 airfoil with angle of attack a = 0°, Re = 0.5 x 10 6 , and M= 0.2. The pressure coefficient C p (x y y) 
at x/L = 0.0, 0.25, 0.5, 0.75, and 1.0, where L is the chord length , and along the surface. The resolution is 383,724 voxels 
and 1,275 surfels for PowerFLOW and 373 x 141 for CFL3D. 
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Fig. 13. NACA-0012 airfoil with angle of attack a = 3°, Re = 0.5 x 10 6 , and M = 0.2. The x-component of velocity 
u(x y y) at x/L = 0.0, 0.25, 0.5, 0.75, and 1.0, where L is the chord length. The resolution is 383,724 voxels and 1,275 surfels 
for PowerFLOW and 373 x 141 for CFL3D. 
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Fig. 14. NACA-0012 airfoil with angle of attack a = 3 0 , Re = 0.5 x 10 6 , and M = 0.2. The y-component of velocity 
v(x, y) at x/L = 0.0, 0.25, 0.5, 0.75, and 1.0, where L is the chord length. The resolution is 383,724 voxels and 1,275 surfels 
for PowerFLOW and 373 x 141 for CFL3D. 










Fig. 15. NACA-0012 airfoil with angle of attack a = 3°, Re = 0.5 x 10 6 , and M= 0.2. The pressure coefficient C p (x y y) 
at x/L = 0.0, 0.25, 0.5, 0.75, and 1.0, where L is the chord length , and along the surface. The resolution is 383,724 voxels 
and 1,275 surfels for PowerFLOW and 373 x 141 for CFL3D. 
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Fig. 16. NACA-0012 airfoil with angle of attack a = 3°, Re = 500 X 10 3 , and M = 0.2. Streamlines of PowerFLOW 
simulation (left) vs. CFL3D simulation (right). 
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Fig. 17. NACA-0012 airfoil with angle of attack a = 3°, Re = 0.5 x 10 9 , and M = 0.2. The x-component of velocity 
u(x y y) at x/L = 0.0, 0.25, 0.5, 0.75, and 1.0, where L is the chord length. The resolution is 383,724 voxels and 1,275 surfels 
for PowerFLOW and 373 x 141 for CFL3D. 
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Fig. 18. NACA-0012 airfoil with angle of attack a = 3°, Re = 0.5 x 10 9 , and M = 0.2. The y-component of velocity 
v(x, y) at x/L = 0.0, 0.25, 0.5, 0.75, and 1.0, where L is the chord length. The resolution is 383,724 voxels and 1,275 surfels 
for PowerFLOW and 373 x 141 for CFL3D. 











Fig. 19. NACA-0012 airfoil with angle of attack a = 3°, Re = 0.5 x 10 9 , and M= 0.2. The pressure coefficient C p (x y y) 
at x/L = 0.0, 0.25, 0.5, 0.75, and 1.0, where L is the chord length , and along the surface. The resolution is 383,724 voxels 
and 1,275 surfels for PowerFLOW and 373 x 141 for CFL3D. 
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